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Abstract — In recent years, the embedding approach for 
switched optimal control problems has been developed in 
a series of papers. However, the embedding approach, 
which advantageously converts the hybrid optimal control 
problem to a classical nonlinear optimization, has not been 
extensively compared to alternative approaches. The goal 
of this paper is thus to compare the embedding approach 
to multi-parametric programming, mixed-integer program- 
ming, gradient-descent based methods, and CPLEX in the 
context of five recently published examples. The five exam- 
ples include a spring-mass system, moving-target tracking 
for a mobile robot, two-tank filling, DC-DC boost converter, 
skid-steered vehicle. A sixth example, an autonomous 
switched 11-region linear system, is used to compare a 
hybrid minimum principle method and traditional numer- 
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ical programming. For a given performance index for each 
case, cost and solution times are presented. It is shown 
that there are numerical advantages of the embedding 
approach: lower performance index cost, generally faster 
solution time, and convergence to a solution when other 
methods fail. In addition, the embedding method requires 
no ad hoc assumptions (such as predetermined mode se- 
quences) or specialized control models. Theoretical advan- 
tages of the embedding approach over the other methods 
are also described: guaranteed existence of a solution 
under mild conditions, convexity of the embedded hybrid 
optimization problem (under the customary conditions on 
the performance index), solvability with traditional tech- 
niques such as sequential quadratic programming avoiding 
the exponential complexity in number of modes/discrete 
variables of mixed-integer programming, applicability to 
affine nonlinear systems, and no need to explicitly label 
autonomous switches with discrete/mode variables. Finally, 
some common misconceptions regarding the embedding 
approach are also addressed including whether it uses an 
average value control model (no), whether it is necessary 
to "tweak" the algorithm to obtain bang-bang solutions 



Febmaiy 27, 2013 



DRAFT 



SUBMITTED TO IEEE TRANSACTIONS ON CONTROL SYSTEMS TECHNOLOGY 2 

(no), whether it requires infinite switching to implement In example (i) we compare the embedding approach 

embedded solution (no), and whether it has real-time ^ Mpp ^jp ^nd CPLEX. For examples (ii) and (iii) 

we compare the embedding method to gradient descent 

based methods in H, ^, and in example (iv) we 

I. Hybrid Optimal Control Embedding ,, , ... , .,, ,, ., , , 

compare the embedding approach with the available 

APPROACH OVERVIEW information in QH and MPP, MIP, and CPLEX 

Hybrid and switched systems have modes of opera- results. In example (v) we show that approach in |[T8], 

tion. Switches can be controlled or autonomous. In much [ 19| is an application of the embedding approach and 

oftheliterature, all such modes are labeled with discrete hence we repHcate their results. Finally, for example 

variables. Solution of these problems is carried out by (vi) we show that the approach in ||9] which reduces 

a number of algorithms including: multi-parametric pro- the autonomous switched problem to a traditional nu- 

gramming (MPP) Q, E), academic E) and commercial merical programming problem favorably compares to the 

(CPLEX) ID mixed-integer programming (MIP) imple- minimum principle solution approach in ^ that utihzes 

mentations, or variants of gradient descent methods H, discrete states to describe the autonomous switches. 
151, 161, 17]. The work shows definitively that the embedding 

In 2005, an embedding approach was developed by method is not only easier to implement, but usually 
Bengea and DeCarlo g) and later extended in Q to faster, achieves universally lower costs, and avoids, in 
a special subclass of autonomous switches. The em- the cases above, the use of simplifying assumptions, 
bedding approach, which converts the switched hy- Specifically, the embedding approach does not require 
brid optimal control problem to a classical nonlinear the use of scheduled switching, the generation of of- 
optimization, has been successfully used in a variety fline maps, and ad hoc assumptions on continuous time 
of appUcations including power management of hybrid control such as a single constant control over extended 
electric vehicles HO) and fuel cell hybrid vehicles lUIl, prediction horizons lUT). For nonlinear affine systems, 
mobile robot slip control IE], and real-time switching the embedding method does not need a piecewise linear- 
control of DC-DC converters O, lE], [IB], lUl]. affine approximation as is done in H], ||l21. Finally, in 

The goal of this paper is to evaluate the convergence the latest implementation of the MPP method 13 the 

time, resulting performance index costs, and require- authors caution against using general nonhnear models; 

ments (e.g., requisite assumptions) of these approaches with the embedding method all that is required is that 

and algorithms in the context of six published examples: the control system be affine. In addition, the embedding 

(i) spring-mass system IT], approach theoretically guarantees the existence of a 

(ii) a mobile robot H, solution assuming appropriate convexity of the integrand 

(iii) two-tank system ||5], of the performance index and an affine control system, 
(iv) a dc-dc boost converter ifTTI . 

(v) skid-steered vehicle lUS], |[T9], and ^- Embedding Approach Description 
(vi) autonomously switched system with 1 1 state space Details of the embedding method have been set forth 

regions 0, Q. in i], M, HI, d, IEQ], IE], El. We first address 
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the relationship of the embedded model to the switched cover the terminal state constraint set with a slighdy 
model, and then switched and embedded optimal control larger open set; this would result in a "new" SOCP 
problems. In the embedding approach to hybrid optimal solution with cost equal to the EOCP solution. The only 
control one forms the embedded model, for example, example in this paper with a terminal constraint set is in 
of a two switched system: i = (1 — v)fo{t,x,uo) + Section Hl-AI a spring mass system, and we construct a 
vfi{t,x,ui). Here w is a control determined to meet switched solution via the embedding approach meeting 
some objective and uq and ui are the continuous time the terminal conditions with a cost lower than MPP 
controls in modes and 1, respectively. In the embedded as published in IH and MIP solutions from the multi- 
model one allows v to take values in the interval [0,1]; in parametric toolbox (MPT) ||2l and CPLEX ||3]. 
the switched model v takes values in the set {0, 1}. Note Thus for the models studied in the examples of this 
that in the switched model, uq and ui cannot be active paper, the EOCP always has a solution and if the 
at the same time and are replaced by a single continuous associated SOCP has a solution, the SOCP solution is 
time control u. One notes that EVERY trajectory of one of the possibly non-unique solutions given by the 
the switched model is a trajectory of the embedded optimal solution space of the EOCP. The SOCP may not 
model. Indeed it has been shown in J8:| that the switched have a solution in the sense that there is no minimum 
trajectories are DENSE in the set of embedded system of the performance index over the set {0, 1} due for 
trajectories. One then defines an integral performance example to constraints: the EOCP solution is the infimum 
index (standard fare) whose integrand is the convex over the set {0, 1} in that case. (Thus the MPP and 
combination of each mode's performance integrand; the MIP approaches are infeasible when there does not exist 
convex combination is analogous to that of the embedded a solution.) Only when there does not exist a SOCP 
model. Each integrand is required to be convex in the solution would one approximate the EOCP solution with 
continuous control. It is important to emphasize that it a switched solution using finite time switching; this is 
has been proved in fS) that the embedded optimal control always possible since the switched system trajectories 
problem (EOCP) always has a solution for any affine are DENSE in the embedded system trajectories. Bengea 
control system assuming the convexity of the integrand and DeCarlo fS] provide a construction for approximat- 
in the continuous-time controls; no similar results exist ing the EOCP solution with an SOCP trajectory to any 
in the hybrid optimal control literature. given precision using the Chattering Lemma; practically 
If the switched optimal control problem (SOCP) has speaking the duty cycle interpretation is often adequate, 
a solution, it is a solution of the EOCP except possibly A brief study of projection methods is set forth in ll23l . 
in the following isolated case: there is a terminal con- All duty cycle interpretations only require finite time 
straint set and an optimal solution reaches the terminal switching. Finally, we point out that when the SOCP 
constraint set at a point on its boundary. In this rare does NOT have a solution, i.e., the performance index 
and easily fixed case, it is possible (but not necessarily does not have a minimum over the class of switched 
true) that the SOCP solution is bounded away from the systems, the costs associated with the embedded solution 
solution to the EOCP; in this case, the SOCP cost is and the projection of that solution via a duty cycle in- 
bounded above the EOCP cost. The "fix" is to simply terpretation onto the feasible set {0, 1}, (PWM solution) 
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are virtually identical |[T3l . ifTSl . |[T6) . original model and simply forms a convex combination 

In general, the SOCP is NOT convex. The EOCP of the vector fields to create a solution space in which 

is convex over rather general conditions set forth in the original problem can be solved. Nowhere is the 

Bengea and DeCarlo fS) which are satisfied by the model averaged in regards to time scales or operating 

formulations in this paper If the SOCP does not have points |[T3], (El, ([15], lfT6l . 

a solution as described previously, then mixed-integer Another question that has arisen is how one guar- 

programming approaches are ill-posed. On the other antees that a bang-bang solution is generated by the 

hand, since the EOCP is convex and there are no integer optimization algorithm? In general, it is not necessary 

variables (even in the presence of autonomous switches) to tweak the algorithm. Non-bang-bang solutions in for 

and since it always has a solution under very reasonable example a 2 mode system require that the Hamiltonians 

conditions, the EOCP can be solved using classical in each mode of operation be equal numerically |J8]; that 

nonlinear programming techniques such as SQP One of situation occurs rarely in general but more often when 

the observations of this paper is that MIP methods fail constraints are imposed on the switching set which must 

evenforsomeof the simple examples tested here; instead remain convex, as in the work of this paper In a 3 

the embedding approach always finds a solution and is mode or greater system such as in [;23i, two or more 

usually faster Hamiltonians can be equal, but if one of the many is 

"larger", then the solution is bang-bang. However, as 

B. Common Concerns about the Embedding Approach ^^^^^^ ^^oyg^ bang-bang solutions are not necessary since 

One common misconception is that the embedding non-bang-bang solutions can always be approximated 

method averages the vector fields of the switched system arbitrarily well with bang-bang solutions, 

similar to the Filippov method in variable structure Another issue is EOCP solutions where v is not in the 

control (VSC). In VSC, Filippov's method is used to set {0, 1}. A value of v that minimizes the performance 

determine a solution to a differential equation whose index in the interior (0, 1) does NOT mean that the 

right hand side is discontinuous on a sliding manifold or SOCP does not have a solution or that infinite switching 

discontinuity surface due to infinitely fast switching in is required to implement the solution. In JS), the value 

the control. Filippov's method takes a convex combina- of v is computed for the example in ||25l which is 

tion of two vector fields (1— a)/++a/^ and chooses the shown to have an infinite number of bang-bang solutions 

variable a to achieve an "average" value consistent with for v = 0.5. In that example, v = 0.5 is shown to 

a tangent plane to the discontinuity surface [241. This is mean that one must spend equal amounts of time in 

not the case in the embedded optimal model where the each mode. Thus a duty cycle interpretation often used 

equivalent of a is a control variable to be chosen so that for a projection method is consistent with KNOWN 

a performance metric is minimized. theoretical properties (|8]. 

Similarly, the boost and buck converter literature con- Further, why is such a solution, v ~ 0.5 , considered 

siders time scale separation and linearization about an singular? Because the Hamiltonians associated with each 

operating point(s) to obtain an "average value model", mode are equal causing a specific function (within a 

On the other hand, the embedding method uses the convex combination of the Hamiltonians) given in JH 
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to be identically zero on a time interval of non-zero approach that additionally finds the mode sequence JS]. 

measure. Equal Hamiltonians and existence of solutions Appendix |A] outlines the MATLAB-based embedding 

were not considered in, e.g., Il26l . Il27l . approach solution algorithm used herein. Space consid- 

A comment similar to others discussed above argues erations preclude detailed descriptions of the optimal 
that the method allows one to use classical sequential control problem solution methods compared to the em- 
quadratic programming, or that the embedding problem bedding approach, however each is briefly presented in 
is nothing more than a classical nonlinear optimization, the example it is first used in. Finally, as set forth in ||9l, 
As pointed out above, that is precisely why one wants traditional numerical methods are applied to an 11 au- 
to use the embedding method, because it completely tonomous mode linear system IS), Q. Passenburg et al. 
avoids the combinatorial complexity of mixed-integer apply a hybrid minimum principle (|6], Q to the problem 
programming. As we will see in the examples to follow, and by formally labeling autonomous modes explicitly 
mixed-integer programming is generally slow and often account for transitions across discontinuity surfaces, 
does not converge. In the following examples, the term "mode" indicates 

In the MPP applications, one does offline computation a dynamical vector field selected with a discrete control 

and then creates look-up tables specific to a model input. However, "mode" has been used in the past to 

and set of objectives. Look-up table approaches cannot indicate both controlled and uncontrolled (autonomous) 

deal with changes in parameters or control objectives, switching of vector fields. Herein, we term the use of 

and hence fail to have the real time capability of the "mode" resulting from an autonomous switch as an 

embedding approach reported, for example, in llT3l . "autonomous-mode" or "a-mode". In ||9l, ET\ . it was 

A conclusion of this work is that for the methods stud- shown that including autonomous switches in the mode 

ied in this paper excluding the embedding approach, not definitions is unnecessary for the embedding approach 

only are there convergence issues for medium horizon and control problems with only autonomous switches 

windows, but when the programs do give an answer, the need no mode designations and are solvable with tradi- 

computation time is possibly orders of magnitude larger tional numerical programming. 

than the embedding method. Further, in each of the examples we use the terms "nu- 
merical optimization cost" and the "simulation cost". By 

II. Hybrid Optimal Control Examples numerical optimization cost we mean the cost computed 

In this section, the embedding approach for hybrid via the numerical optimization program using collocation 

optimal control (HOC) is applied to a spring-mass hybrid and trapezoidal numerical integration of the PI; by simu- 

system HI, a switched-mode mobile robot lIU, a two- lated cost we mean the cost obtained by numerically in- 

tank hybrid system (5], a DC-DC boost converter ifTTI . tegrating the system ODE's using the piecewise constant 

and a skid-steered vehicle ifTSl . lfT9l . The results from continuous controls from the numerical optimization, 

the embedding approach are compared to those from In all of the examples, the EOCP is solved using MAT- 

MPP m, MIP/MPT (MIP using the MPT), MIP using LAB 's>i;nco« function following the general procedure 

CPLEX, a method that computes switching times for a outlined in Appendix [A] Using a numerically superior 

pre-determined mode sequence |J4] as well as a follow-on solver would only improve the convergence rate and 
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solution times of the results reported for the EOCP. All 
the code used to generate the results in this paper can 
be accessed at 



A. Spring-Mass Hybrid System ^7]/ 

Example 14.2 in HI introduces the spring-mass hybrid 
optimal control problem. A mass is connected to ground 
with a spring in series with a damper that represents 
viscous friction. The spring has affine characteristics and 
the viscous friction coefficient can be changed from one 
value 61 to a different value &2 instantaneously with a 
binary input. The continuous-time spring-mass system 
dynamics are 

ilit) =X2{t) 



Mx2{t) = - k{xi{t)) - b(u2{t))x2{t) + Ui{t) 



(1) 



where the the spring coefficient k{xi{t)) and viscous 
friction coefficient b{u2{t)) are 



fcixi(i) +di, xi(t) < XjTL 
Hxiit))={ (2) 

k2Xl{t) +d2, Xi{t) > X,n 



HMt)) 



61, U2{t) = 1 



(3) 



62, U2{t)^0 

with xi and X2 the mass position and velocity, respec- 
tively; the system has two modes and two a-modes. 
System parameters are given as M = l,bi = 1, 62 = 50, 
fci = 1, fc2 = 3, di = 1, d2 = 7.5, and x„i = 1- The 
spring-mass example is solved using MPP (the original 
solution approach in U]), mixed-integer programming 
with both MIP/MPT and CPLEX, and the embedding 
approach. MPP IT], ||2| solves the optimization problem 
offline for a full range of initial state values to obtain 
explicit state-feedback controllers. The optimal state- 
feedback control for a linear piecewise affine (PWA) 
control model (as used with MPP here) is a PWA 



function of the initial state with controller parameters 
defined on polyhedron partitions of the feasible state 
space, i.e., the controller parameters are taken from a 
look-up table with initial state input. 

To implement the MPP algorithm in |[T|, four "modes" 
were defined depending on the discrete input U2 and 
spring coefficient of (|2]i; the "modes" are a mixture of 
both modes in the sense used here and a-modes. The 
dynamics in each "mode" are discretized with a 0.5 time 
unit sampling resulting in a discrete-time PWA system: 

x{k + 1) 

f 

Aix{k) + Biui{k) + fi, .Ti(fc) <l,U2(fc) <0.5 

A2x{k) + B2Uiik) + f2, .Ti(fc) > l,M2(fc) < 0.5 

A3x{k) + B^uiik) + h, xi(fc) < l,U2{k) > 0.5 
A4x{k) +B4Ui{k) + /4, xi(fc) > l,U2(fc) > 0.5 

(4) 

where Ai, Bi, and fi are listed in [Qj. The MPP control 

objective is to minimize 

JraixiO),Uo,N)^x{NfPx{N) 

N-1 

+ Y, x{kfQx{k) + u{kfRu{k) 



fc=0 



(5) 



over Uo = [uiOY , . . . , u{N - l^f (u{k) = 
[ui{k),U2{k)]'^) and subject to (|4]i, xi,X2 € [—5,5], 
Ml e [—10,10], and a terminal state constraint set 
Xf e [-0.01,0.01] X [-0.01,0.01]. The performance 
index parameters are A^ = 3, P = Q = I2, R = [ V il 
and the initial state is a;(0) = [3, 4]^. The MPP problem 
is solved using MPT routines in MATLAB. Furthermore, 
the MPP problem formulation is also solved with both 
MIP/MPT and CPLEX. The MIP/MPT problem is also 
solved with MPT routines in MATLAB. For the CPLEX 
problem, CPLEX's convex quadratic performance index 
optimizer, subject to linear constraints, is called from 
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MATLAB. 

Next, the embedding approach problem is set forth. 
First, the continuous-time embedded model is formed 
using two modes (recall the embedding approach does 
not define separate modes based on autonomous switch- 
ing 191, i.e., have a-modes). 



i{t) ={i - v{t)) 



■k{ii) - b2X2{t) + u\{t) 

i2{t) 
k{5:i) - bii2{t) + u\{t) 



(6) 



where (~) is an embedded system value, v G [0,1] is 
the embedded mode switch value and u\ and u\ are 
the continuous controls in each mode. Notice that v = 
0/v = 1 in (|6]l corresponds to U2 = 0/u2 = 1 in ([TJ. The 
embedded problem is 



min j£;(a;(<g),Wi,'Ui,w,ig,i?) 



where [^O'^^l ^^ ^^^ prediction interval and 



(7) 



(8) 



Je 






{i{tfQE{t)i{t) 



+ {i-m) 



u°{t) 



+i{t) u\{t) 1 RE{t) 



RE{t) 

ulit) 
1 



SO(i) 



dt, 



(9) 



the embedded context as 

i{k + 1) 

' (1 ~ v{k)) [Aii(fc) + Bi«0(fc) + /i] 
+i{k) [A^iik) + Bsulik) + h] , 

{l-i,(k))[A2i{k)+B2U°,{k)+f2] 

+vik) [Aii{k) + Biulik) + U] , 



Xi 



(fc)<l 



Xl(fc) > 1 

(10) 
and the trapezoidal numerical integration representation 
of© is 

N-l 

Je ^i{NfPS:{N) + ^ {i{kfQi{k) 



(11) 







k 


=0 








(i-Hk)) 


u?(/s) 


R 


u",{k) 



v{k) 


ul{k) 1 


R 


u\(k) 
1 


1 



with the minimization subject to ^, the previously 
defined bounds on the states and continuous control 
input, and terminal state constraint set. To solve the 
embedded optimal control problem numerically, the dis- 
cretized system representation in (@) is reformulated in 



where QE{t — to) and REit — t^) were chosen such that 
( fTTT i is equivalent to ^. If v{k) e (0,1), then mode 
and control projection are required. Here, the projection 
method of lITTI . ll23l is applied where the projected 
mode, V, is zero if w < 0.5 and one otherwise. The 
projected control ui is equal to (1 — ?})u5 if w = and 
im\ if {> = 1. 

Figs. [Hand 12] shows the state trajectories and evolution 
of the continuous and discrete controls for the embed- 
ding approach and MPP; MIP/MPT and CPLEX results 
are nearly identical to the MPP results and not shown. 
The embedding approach results in 3 mode projections 
over 25 partitions. The projections occurred on the 
intervals [0.5, 1], [1, 1.5] and [12, 12.5]. In all cases, the 
embedded mode value was less than approximately 0.3, 
which implies that the mode value was projected to 0. 
Table U lists the performance index costs, ,/, and control 
solution times, Af, for the simulation associated with 
A^ = 3 plus simulations with A^ = 4 with ts = 0.375 
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Fig. 1. Spring-mass example state responses: ( — ) EOC with projec- 
tion, (•) MPP (MIP/MPT and CPLEX liave similar results). 
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Fig. 2. Spring-mass example continuous and discrete controls: ( — ) 
EOC with projection, ( ) EOC without projection, (•) MPP (MIP/MPT 
and CPLEX have similar results). 



and iV = 5 with ts = 0.3. There are 10 mode 
projections over 33 partitions needed for A^ = 4 and 
7 over 42 partitions for A^ = 5. We conclude that for 
the same performance metrics and the same example, 
the embedding approach achieved lower costs in all 
trials. Furthermore, for a model predictive control (MPC) 
window of 5 partitions, the embedding approach was 
143 times faster than the MIP/MPT approach and 22 



TABLE I 

Spring-mass example EOC (with and without proiection), 

MPP, MIP/MPT, CPLEX PERFORMANCE INDEX COST, J, AND 
SOLUTION TIMES, Af, FOR THE SIMULATION. 





N 


= 3 


/V 


= 4 


N 


= 5 


Metric 
















ts = 


= 0.5 


ts = 


0.375 


ts = 


= 0.3 



Je, proj. 


Aim 


62.38 


85.65 


Je, no proj. 


47.50 


59.70 


71.67 


J MPP 


53.90 


72.84 


89.35 


Jmip/mpt 


53.81 


72.78 


89.35 


JCPLEX 


53.89 


72.84 


87.73 


At^j, proj. 


3.86s 


9.42s 


9.99s 


At^, no proj. 


3.81s 


7.78s 


8.75s 


'-^^MPP 


9.31s 


31.20s 


222.55s 


Af 

^''MIP/MPT 


33.73s 


72.78s 


1434.7s 


^icLEX 


1.00s 


2.12s 


3.20s 



times faster than the MPP approach. Although the MPP 
approach generates a look-up table, if any parameters 
change or there are modeling issues, the entire table 
must be recomputed. The embedding approach was on 
average about 3.5 times slower than using CPLEX. 
CPLEX is not a general purpose minimization program 
like finincon that we use to solve the EOCP. CPLEX is 
specifically tailored to minimize a quadratic cost function 
subject to linear equality constraints (as is the spring- 
mass problem). 

B. Mobile Robot Hybrid System JSJj 

Wardi et al. ID consider the problem of a mobile robot 
tracking a moving target while avoiding obstacles. The 
robot has three operating modes: go to goal (G2G), avoid 
obstacle 1 (Avoidl), and avoid obstacle 2 (Avoid2). The 
robot's dynamic equation of motion is 



m 



FcOs(x3(i)) 

Vsm{x3{t)) 

U-X3{t) 



(12) 



Febmary 27, 2013 



SUBMITTED TO IEEE TRANSACTIONS ON CONTROL SYSTEMS TECHNOLOGY 



where V is 



V ^ 



V, 



\xr - xgW > r 



(13) 



^\\xr-xg\\, \\xb^ -XgW < r 
and u is defined for the three modes as 

UG2G = tan ' 



'^Avoidl 



(14) 



(15) 



(/>$ — 7r/2, 0* — 2:3 > 

0$ + 7r/2, (f)^ — X3 < 

(f)xj, — 7r/2, 0* — 2:3 > 

UAvoid2 = •\ (16) 

(f'^i + 7r/2, 0* — 2:3 < 



and 0$ and (/>$ are 

0$ = tan 
0$ = tan 



1 fx^^2 -'XR^2{t) 



X<i>,l -XR^i{t) 
-1 (x^,2 -~XR^2{t) 



(17) 
(18) 



^x^.i - 2:fl,i(t) 

where x ~ [xr^X'^Y' such that xr ~ [xb,.i^xr^2V is 
the global position (two coordinates) of the mobile robot 
and xs G [0, 2tt) is the robot heading angle; xg is the 
coordinate pair of the goal; V^ is a robot constant speed 
equal to 1; r is 0.5; a;$ = [0,4]^ is the coordinates 
of obstacle 1; and x^, ~ [6,6]^ is the coordinates of 
obstacle 2. 

The control problem given is to select the G2G, 
Avoid 1, or Avoid2 that minimizes 



J{x{t),XG{t),t) 



p\\xR{t)-XG{t)f 



+ ai exp 
+a2 exp 



\\xR{t)-x^f 

Pi 
WxRJt) -2:vi,|p 

/32 



dt 
(19) 



subject to the system equations from the current time t to 
the end of the simulation tf. The cost function constants 
are p = 0.1, ai = 03 = 500, and /3i == /32 == 0.8. 
The goal coordinates are not known a priori, rather an 



estimate, xg at s > Hs used in the cost function: 

xg{s, t, XG{t)) := XG{t) + XG{t){s - t) (20) 

where the derivative of XG{t) is approximated with 
[xcit) — XG{t — ts))/ts. The problem is defined as 
having to = 0' ^/ = 18, and a sample time, <s, of 0.1. 
Wardi et al. |J4| solve the control problem by finding the 
switching times for a pre-determined mode sequence. 
During a simulation, a gradient-descent technique is 
used to regularly update the mode switching times for 
upcoming modes left in the sequence. Their results are 
reported further on. 

MPP, MIP/MPT, and CPLEX are not suitable for 
problems having nonlinear dynamics and constraints 
and/or Pis not linear or not quadratic as in this example 
of mobile robot tracking and avoidance. The embedding 
approach has no such limitations. Hence only the embed- 
ding approach is used for comparison with the published 
results in ID. 

The three modes carry over to the embedded model 
and are controlled by wo, wi, and V2'- 

V^cos(£3(t)) V^cos(£3(i)) 

x{t)^vo{t) T/sin(x3(t)) +ii{t) V sm(ii{t)) 

UG2G - Xs,{t) UAvoidl - Xs,{t) 

1/C0s(x3(t)) 

+W2(t) ysin(i3W) 

UAvoid2 — X'i{t) 

(21) 
where u,; G [0, 1] and vq + ^1 + ^2 = 1- The em- 
bedded control problem is to also minimize (fT9] l over 
the Vi subject to ^T\\ and the summation constraint on 
Vi. MATLAB's fmincon function is used to solve the 
problem after converting the continuous-time equations 
to discrete-time equality constraints via collocation ||9l. 
The cost function is approximated with trapezoidal nu- 
merical integration (see Appendix |A] for more detail). 
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Fig. 3. Mobile robot trajectory over the simulation: ( — ) EOC with 
projection, (•) EOC no projection, (■) Wardi et al. off-line solution (4], 
(■) target path, (o) obstacle 1, (D) obstacle 2. 
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Fig. 4. Mobile robot embedded optimal control mode selection: ( — ) 
projected, (•) embedded solution. 



A continuation approach 11291 is used to construct an 
initial guess; furthermore, the continuation approach 
itself is initialized using Appendix [Aj-Step |4l When 
needed, mode projection is accomplished by selecting 
the embedded mode switch with the largest value as was 
done in 1231 . 

Fig. [3] shows the robot trajectories obtained with the 
embedding approach and using the switching times and 



mode selection given in ID for the off-line solution. 
Fig. m shows computed values for the embedding ap- 
proach and final projected mode switch values. Mode 
projection is required approximately 87% of the time. 
In Wardi et al., the mode sequence is predefined as 
{G2G,Avoidl,G2G,Avoid2,G2G} and only the switch- 
ing times are to be found. The embedding approach 
makes no such assumptions on the mode selection. The 
cost obtained with the Wardi et al. offline solution is 
22.84 while the embedding approach here results in a 
simulation cost of 18.83, a reduction of nearly 18%. As a 
follow-on test, an MPC approach was implemented with 
a prediction horizon of 1 time unit (10 partitions of 0.1 
time units). The MPC optimization and implementation 
resulted in a simulation cost of 20.56, which is still about 
10% less than the offline solution cost of 22.84, and 154 
mode projections (out of 180 partitions). 

C. Two-Tank Hybrid System ^5^ 

The two-tank hybrid system in |]5] consists of the 
draining of one tank into another with the flow rate 
into the first tank regulated by a valve with two possible 
values, ly = 1 and ly ~ 2. The dynamics are given as 



i(f) 



z/(t) - y^T(i) 



(22) 



where x = [xi,a;2]"^ with xi the level of tank 1 and X2 
the level of tank 2. The initial conditions are .ti(0) = 
X2(0) — 2 and the states satisfy a;i,a;2 G [li4]. The 
control problem is to select the i'{t) that minimizes 

/•*/ 
J{x{0),tf)=2 {X2{t)-3fdt (23) 

Jo 

with tf = 20 and subject to ( l22l i and the state bounds. 
The sample interval is tg = 0.01 and the continuous- 
time dynamics are discretized with the forward-Euler 
method. Wardi and Egerstedt ||5] use a gradient-descent 
technique to solve the control problem like that used for 
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the mobile robot described previously except the need 
to predefine the mode sequence is removed. The mode 
sequence and switch times are calculated to minimize the 
cost function over the entire simulation time, i.e, it is a 
global minimum solution approach. The approach relies 
on an insertion gradient that indicates whether a change 
in the mode sequence and/or switch times decreases a 
cost function. The insertion gradient is iteratively driven 
toward zero over mode sequences and switch times 
selected according to the solution algorithm. 

The two-tank problem is solved using the embedding 
approach. The MPP, MIP/MPT, or CPLEX methods ai-e 
not applied to the two-tank problem because of their 
inability to handle models with square-root terms. The 
embedded model representation of the two mode two- 
tank dynamics is 




10 12 14 16 18 



20 



4 
3 

2 

I 

1 

0.5 





10 12 14 16 18 20 



m=ii-m) 



i{t) 



l-^iAi) 



\/xiit) - \/x2it) 

2- V^TW 



(24) 



V^i(i) - \/x2 (i) 

where v G [0, 1]. Notice that the embedded problem does 
not include continuous control inputs, only the mode 
needs to be selected. The embedded control problem is 
to minimize ( |23] ) over v subject to (l24l i and the bounds 
on the states listed previously. Appendix |A] outlines 
how MATLAB's finincon function is used to solve the 
problem after converting the continuous-time equations 
to discrete-time equality constraints via the forward- 
Euler method and approximating the cost function with 
trapezoidal numerical integration. Mode projection is of 
the duty cycle type here. For (1 — v)ts^E, mode with 
i^ = 1 is applied and then for the remainder of the sample 
interval mode 1 with z^ = 2 is used; ts^s is the embedded 
problem sample interval. 

Fig. |5]shows the embedding approach performance for 



2 4 6 8 10 12 14 16 18 20 
Time 

Fig. 5. Two-tank hybrid system example states and embedded 
mode switch values: ( — ) EOC with duty cycle mode projection, (•) 
embedded solution. 



ts.E ~ 0.25. The X2 state reaches the desired value of 3 
about 6 time units into the simulation and does not vary 
more than 1.5% thereafter In contrast, in [0 the state X2 
shows variation up to about 7% from 6 time units onward 
and both states show signs of instability as t approaches 
tf. The reported cost is 4.78 (recall ts ~ 0.01). On 
the other hand, for a time step of 0.25, the embedding 
approach gives a lower (simulated plant) cost of 4.74 
and shows no signs of instability. 

D. DC-DC Boost Converter Hybrid System M7'i 

The DC-DC boost converter model is a switched 
lumped parameter circuit described in IITtI . The boost 
converter continuous-time dynamics on the time interval 

[fcis, (fc + l)is] are given as 



x{t) 



Fix{t) + fiVs, kts<t<{k + d{k))ts 





l^ 


2x{t)+f2Vs, {k 


+ d{k))ts <t <{k + l)t 


ith 


(25) 


i^i 


^H, 


n 1 


H{' (26) 




L ^ C(r„+r,)J 
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Hi 


= 
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To 






(28) 




To+ra. 






1 f^. ^ '■or-c ^ 


ro 




F2^H2 




'' ' To 
■'■0 


fr.y 


1 


i?2"' (29) 




^{ro+rc) 


C(ro+re)J 




f2=H2 


i 0; 


T 




(30) 


H2 


= 


1 

_ro+rc To 




To 




(31) 




+ra_ 





where x = [ii,Vo\'^, ii is the inductor current, Vo is the 
output vohage, Vg is the supply voltage, and d G [0, 0.95] 
is the duty cycle in the interval. The model parameters 
ai-e C = 100 fiF, L = 2 mH, r^ = 0.1 Q, n = 0.5 VL, 
To = 200 fl, and the maximum inductor current is 2.5 A; 
the switching frequency is 20 kHz (the sample interval, 
ts, is 50 /is). 

For MPP control in ifTTl . ( l25T l is not used in the actual 
optimization. Rather, linearized, discrete-time models 
developed from (l25l l corresponding to one of three 
regions of the duty cycle interval, Dq = [0.0.45], 
Di = [0.45,0.6], and D2 = [0.6,0.95], are used: 

x'{k + 1) = A,„,,x'(fc) + B,„,,d(fc) + /„,,, (32) 

where x'{t) = [i'i,iQ'^ = [ii/vs,Vo/vs]'^ (the original 
state is divided by Vg to avoid the need to generate new 
Unear models if Vg changes) and i = 0, 1, 2 designates 
the duty cycle region of the model. AppendixlBJdescribes 
the process of generating the linear model in IITtI and 
lists the A„i_i, B„i,i, and fm.i found. 

The MPP control performance index given in |[T7) for 
a horizon window of length A^ is 

N-l 

J{x{k),Dk,N) = Y, WQKik + l\k)^ v'o.refimi 

k=0 

+ \\Rid{k + l\k) - dik + I - l\k))\\l 

(33) 



where Dk = [d{k ~ 1), . . . ,d{k + N - 1)], Q and R are 
penalty weight scalars, v'^ re/(^) ~ "^^o.ref/vs is a given 
reference output voltage (assumed constant over A^), and 
p designates the norm. The problem parameters given 
are x{0) = [0, v^f, N = 18, Vo,ref = 50 V, w, = 25 V, 
Q = 10, i? = 1, p = 1, and q = I. The application of the 
duty cycle to the plant is given as being delayed by one 
sample time interval. In ifTTl . the duty cycle is fixed to 
the same value over the A^ partition prediction horizon, 
reducing the optimization to a simple one-dimensional 
search rather than an A^-dimensional search as would be 
consistent with a normal MPC framework. 

The MPP, MIP/MPT, and embedding approach are 
compared using a 1-norm PI as in the original problem; 
CPLEX and the embedding approach are compared using 
a squared 2-norm (quadratic) PI because CPLEX is not 
able to calculate a p-norm PL Also, the limiting assump- 
tion of a fixed duty cycle over the prediction horizon is 
removed. Furthermore, the embedding approach theory 
is not yet developed for p-norm weights on the embedded 
switch values, i.e., the duty cycle values as seen later, 
thus R is set to zero in the subsequent HOCP solution 
approaches. 

The MPP MIP/MPT, CPLEX and the embedding ap- 
proach algorithm specific boost converter HOCP formu- 
lations are now set forth. To apply MPP and MIP/MPT 
control using the MPT, x'{k) is augmented with d(fc— 1) 
and Wg^gf states to form a new state vector x"{k) = 
[i'i{k) , v'g{k) , d{k — 1), w^^gf(fc)]^; the control input is 
Ad{k). The equivalent performance index to (|33) is 

N-l 

J{x"{k),ADk,N) = Y. \\Q'^"{k+l\k)\\i + \\RADk\\i 

1=0 

(34) 

where ADk = [Ad(fc), . . . , Ad{k + N - 1)]^, R ^ 0, 

and Q'(2,2) = 10 and g'(2,4) = -10. 

CPLEX is the next solution approach to the boost 
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converter HOCP presented. CPLEX is not able to solve 
problems with p-norm performance indices therefore the 
quadratic PI is ( |33] | with p = 2 and q = 2. Also, the 
boost converter dynamics are modeled in CPLEX using 
( |32] |, developed for the MPP and MIP/MPT control, 
because CPLEX accepts only linear equality constraints. 
The embedding approach is now considered for the 
boost converter HCOP . Unlike for MPP and MIP 
(MIP/MPT and CPLEX), the ti'ue boost converter dy- 
namics, dZSl ), are used in the optimization; there is no 
need for ad hoc pseudo-model generation. The boost 
converter dynamics in the embedded model formulation 



TABLE II 

DC-DC BOOST CONVERTER EXAMPLE EOC (DUTY CYCLE 

INTERPRETATION), MPP, AND MIP/MPT 1-NORM PERFORMANCE 

INDEX COST J, AND SOLUTION TIMES, Af, FOR THE SIMULATION. 

(X IS NO DATA.) 



Metric 




N = A 


Af = 6 


Af = 12 


Je 




917.33 


818.45 


808.94 


J MPP 




5651.93 


X 


X 


Jmip/mpt 




970.09 


954.95 


X 


Af^ (s) 




85.17 


804.50 


3789.53 


^*A/PP (S) 




2630.98 


X 


X 


'^''MIP/MPT 


(s) 


103.87 


1145.07 


X 



(35) 



x{t) ^{l - i{t)){F23:{t) + f2Vs) 
+ vit){Fii{t) + fivs) 
The embedded optimal control problem is to minimize 

JE{xitl),V,tltP) ^ I \\QE{io{t)-Vo,ref{t)Wpdt 

(36) 
over V <S [0, 1] subject to ( |35] |. Qe, P, and q are 
chosen such that the trapezoidal numerical integration 
approximation of (|36] | is consistent with the PI of the 
MPP, MIP/MPT, or CPLEX approach being compared. 
The EOCP is solved using MATLAB's^mco« function 
after converting the continuous-time dynamic equations 
to discrete-time equality constraints via collocation and 
approximating the cost function with trapezoidal numeri- 
cal integration (see Appendix lAl for the solution method 
outline). Mode projection of an optimal v consists of 
using it directly as a duty cycle, d over a partition. 

The comparisons of the boost converter HOCP solu- 
tion methods is carried out using iV = 4, A^ = 6, and 
A^ ~ 12 and reference output voltage of 



TABLE III 

DC-DC BOOST CONVERTER EXAMPLE EOC (DUTY CYCLE 

INTERPRETATION) AND CPLEX SQUARED 2-NORM PERFORMANCE 

INDEX COST WITH (36), J, AND SOLUTION TIMES, Af, FOR THE 

SIMULATION. 



Metric 



Af : 



Af = 6 



Af = 12 



Je 




12593.16 


9969.59 


9291.08 


JCPLEX 




156547.91 


11512.15 


10297.26 


At^j (s) 




50.79 


81.82 


267.17 


^^CLEX 


(s) 


96.06 


129.56 


31605.55 



'Jo,ref{t) 



35 y, 25 ms < t < 35 ms 



50 V, otherwise. 



(37) 



Longer prediction horizons, like A^ = 18, are not 
considered since it has been shown that longer prediction 
horizons do not provide significantly better control ll30l . 
Table M lists the MPP, MIP/MPT, and embedding ap- 
proach 1-norm PI costs and total control solution times. 
In all tests, the embedding approach provided the least 
cost and fastest solution time when a comparison could 
be made to MPP and MIP/MPT. In several instances, 
MPP and MIP/MPT solutions are not listed because 
control solutions were not generated after 40 hours and 
the test was stopped. Table |III] shows the CPLEX and 
embedding approach quadratic PI costs and total control 
solution times. The embedding approach costs are lower 
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Fig. 6. Boost converter hybrid system states and duty cycle from 
embedding solution approach with A^ = 12: ( — ) EOC with duty cycle 
mode projection, (•) CPLEX, ( ) I'o^ref- 



than the CPLEX costs in each test. Also, as TV is 
increased the cost decreases for both solution methods. 
The embedding approach solution time for iV = 4, 
N = 6, and N = 12 are about 2, 1.6, and 120 times 
faster than the CPLEX times, respectively. The fact that 
the CPLEX solution time dramatically increases over 
the embedding approach time when A^ = 12 illustrates 
the "curse of dimensionality" associated with mixed- 
integer programming. Boost converter control using the 
embedding approach implemented in a dedicated micro- 
controller has been solved as fast as 50 ^s |[T3l . 

Fig. |6] shows embedding and CPLEX approach results 
for A^ = 12. The ( l37b reference output voltage appears 
to be better tracked with the embedding approach during 
the step change in Vo,ref to 35 V on [25, 35) ms. The 
embedding approach is better able to track the 50 V 
reference on [10,25] ms (after the initial voltage rise 
and before the step change in vq) than CPLEX; the 
average Vo on [10, 25] ms from the embedding approach 
is 50.15 V while that from CPLEX is 49.69 V, the 
error with the embedding approach is about half that 



of CPLEX. 

One final point, good control performance requires an 
accurate value of the possibly varying load resistance 
and good estimates of the resistance have been achieved 
using a nonlinear resistance error system [1T61 . MPP, 
MIP/MPT, nor CPLEX has the capability to incorporate 
a nonlinear estimator directly into the control problem 
like the embedding approach. 



E. Skid-Steered Vehicle Hybrid System I118]I , 4791/ 

A skid-steered vehicle (SSV) movement hybrid prob- 
lem is set forth by Caldwell and Murphey in ifTSl , 
||T9l . They identify four modes of operation during SSV 
movement depending on the lateral sticking or sliding 
of the wheels; mode 1, front and rear wheels sticking 
laterally; mode 2, front wheels sticking laterally and back 
wheels skidding laterally; mode 3, front wheels skidding 
laterally and back wheels sticking laterally; and mode 4, 
front and rear wheels skidding laterally. The equations 
of motion for modes 1 and 4 are given in il9l'l ; mode 2 
and 3 dynamics are fisted in Appendix |C] due to their 
length. 



Xit)- 
Y{t) = 



(F1+F2+F3+F4) 


cose(t) 


M 




(F1+F2+F3+F4) 


3ine(t) 


M 





- ciX{t) 
ciY{t) 



(38) 



'The sign on the 9 terni in inode 4 is negative here rather than 
positive as in 1191 and was changed after con'espondence with the 
authors. 
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X{t) = (f-l+J-2+J-3 + F4)cose(t) _ ^^^(^^ 

+fikg sin e{t) [-X{t) sin e{t) + Y{t) cos 0{t)] 



switches. The SSV HOCP is to minimize the following 
cost EH: 



Y{t) = i^^+^2+^^+r,,s,n,(T) _ ^^yf^^^ 



{Fi+F2+F3+Fi)sine(t) 
M 

Hk9 cos 9{t) [~X{t) sin e{t) + Y{t) cos 9{t)] 



m 



_ b{Fi-F2-F3+Fi)-a^tJ.KMge{t) 
~ J 



J{V) =0.5(5;(t/) - Xref{tf)fP{i{tf) ~ Xrefitf)) 

0.5{x{t) ~ Xref(t)fQ{x{t) - Xref{t))dt 

(42) 



to 



with 



M =mb + Am^ 



Am,, 



J =4m^(a2 + b^) + -j^iwl + 3wl) 



(39) 



(40) 



(41) 



subject to a relaxed switched system representation of 
the dynamics, i.e., an embedded representation: 



i=l 


(43) 




(44) 



where X, Y, and 9 are position and orientation co- 
ordinates in the global coordinate frame; F1/F2 is the 
right/left rear wheel torque; F4/F3 is the right/left front 
wheel torque; a are mode specific damping coefficients; 
fiK is the coefficient of kinetic friction; mi, is the body 
mass; m^u is a wheel mass; BilBu, is the SSV body 
length/width; alb is half the distance from wheel center 
to wheel center in the vehicle length/width direction; 
uitu is half the wheel width; and Wr is the wheel 
radius. Model parameters from llTSl . |1T91^ are ci = 0.7, 
C2,C3,C4 =: 1.2, ^iK = 0.8, mb ^ 70 kg, m^, = 2.5 kg, 
Bi = 0.8 m, B.^ = 0.6 m, a ^ 0.16 m, b ^ 0.28 m, and 
Ww,Wr = 0.14 m. 

Caldwell and Murphey use the SSV mode-dependent 
model to first generate the vehicle response given wheel 
torques (shown in Figure 3 in |fT9l ) and mode values 
over a 15 s period. The mode sequence is (1,4, 1) with 
switchings occurring at 5 s and lis. Next, the simulation 
output is used as reference values to be tracked in a 
HOCP in an effort to try to recover the simulation mode 

-Values for C2, C3, w^, and Wr were obtained through personal 
correspondence with Caldwell and Muiphey. 



where i = [X, X, Y, Y, 9, 9f, Q = diag(l, 1, 1, 1, 1, 1), 
P = diag(l, 10, 1, 10, 1, 10), to = s, and tf = 15 s. 
Finally, the active mode is chosen as the mode associated 
with the greatest value of Vi. This is the mode projection 
method in Meyer et al. 1231 found to give the least cost 
in an empirical study of projection methods applied to 
an example hybrid optimal control problem solved with 
embedding. Caldwell and Murphey lfT9l showed that the 
embedding approach and mode projection can recover 
switched system modes for the SSV described here. 
Their approach was presented as a heuristic development 
independent of the work in f8|. In contrast, E provides 
a theoretically rigorous justification of the embedding 
approach and shows that it essentially always leads to a 
solution of the SOCP. 

The embedded problem is solved here and the mode 
projection method used in Caldwell and Murphey is 
appUed. Again, due to the nonlinear nature of the prob- 
lem, MPP, MIP/MPT, and CPLEX methods cannot be 
utilized. For the embedded problem, the 15 s simu- 
lation time is divided into 60 equal length partitions. 
The HOCP embedding approach solution follows Ap- 
pendix |A] where the continuous-time cost integration is 
represented with trapezoidal numerical integration and 
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Fig. [U Each region also has an associated PI integrand: 

Fg{x, u) = 0.5x^{t)Sqx{t) + 0.5u^{t)Rqu{t) (46) 



where Sq,Rq G 



^^^ are constant penahy weight 
matrices that depend on q. Given a system tra- 



5 10 

Time (s) 



15 



Fig. 7. Skid-steered vehicle embedded solution and projected modes: 
( — ) embedded, (•) projected. 



embedded continuous-time dynamics are transformed 
into equality constraints using collocation. The original 
mode sequence and switching times used to generate the 
reference state trajectories is recovered using embedding 
and projection as shown in Fig. |2l mode projection was 
required for 10% of the partitions. Also, the embedded 
solution cost is 2.1-10^'', meaning the embedded solu- 
tion tracked the reference trajectories well. 

F. 11 State-Space Region Autonomous Switch Exam- 
ple 1^ 

Passenberg et al. ||6l, Q proposed a version of a 
hybrid system minimum principle and applied it to an 
1 1 a-mode linear system. In the example, a state space, 
M?, is divided into 1 1 polygonal regions, each associated 
with a distinct linear state dynamic. Thus, a discrete state 
q{t) G {1,2,..., 11} identifies a continuous-time state 
dynamic: 

x{t) = Aqx{t) + Bgu(t) (45) 

where x eR'^, u e Uq CR^, Aq eM'^^^ constant, and 
Bq e M^^^ constant. The 11 regions are illustrated in 



jectory {x{t),u{t)), the associated discrete state se- 
quence {qi, . . . jQk} and the switching time sequence 
{to,...,iif} {to and tK are the initial and final time, 
respectively), such that x{t) belongs to region qj if 
tj-i < t < tj, the PI is then: 

^ ptk 
J{x{t),u{t))=Y. ^9^-^^- (47) 

The optimal control problem is to drive the states 
from x(fo) = [—8,-8]^ at tg = to the origin at 
tK = 2 such that the PI in (07]) is minimized subject 
to ( |45] | and constraints on u due to Uq (which is not 
the same for all values of q). The Aq, Bq, Sq, Rq, 
and Uq are listed in Appendix |D] along with the state- 
space partition boundary equations. In this 1 1 state-space 
regions problem, all the switches are autonomous. At 
the switching surfaces both the vector fields and the 
performance index are nonsmooth, which is problematic 
for traditional optimization methods such as SQP 

The optimization algorithm in Q can be summarized 
as follows: 

1) Choose a feasible sequence of discrete states (adja- 
cent regions) and an appropriate number of switch- 
ing points (switching time and the value of the state 
at the switch) on the switching manifolds. 

2) Find the optimal value for the continuous control 
u{t) by solving a set of traditional optimization 
problems for which the switching points provide 
the fixed initial and final states and times. 

3) Determine the gradient of the PI at the (au- 
tonomous) switching points. 
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4) Based on the gradient computed in Step [3] and the 
local geometry of the switching manifolds around 
the switching point, determine the next sequence of 
the discrete states (regions) and switching points. 

5) Stop if the next switching points are within a certain 
tolerance of the current ones otherwise return to 
SteplH 

In our view. Step |4] is the critical step and the main 
insight of Passenberg et al. In particular, by separating 
the computation of the continuous controls from the 
computation of the switching points, the nonsmoothness 
at the switching points is completely avoided. Passen- 
berg et al. |[6| reported a solution cost of 14.69 for 
the example with 172 nonuniform time partitions; the 
solution trajectory is shown in Fig. QD- 

In our previous work |]9], ET\ we showed that prob- 
lems involving a-modes, i.e., autonomous switches, can 
be solved using traditional mathematical programming 
methods (in particular, SQP) without any special consid- 
eration of the switches and the nonsmoothness associated 
with them. However, if applied to the example above, the 
direct appUcation of SQP results in a slightly suboptimal 
solution. 

To evaluate our approach, we first applied traditional 
numerical programming as suggested in ||9] with 172 
uniform length partitions; for better numerical solution 
behavior, Bq and Rq were scaled so that the continuous 
controls, u, lie in [—1,1] x [—1,1] for each q. Since 
all switches are autonomous, no embedding in the vein 
of II2TII is necessary. Starting from an initial guess 
obtained using continuation 1291 (which is initialized 
using Step |4] in Appendix |A|, we obtained a numerical 
optimization cost of 14.71 and simulated plant cost with 

^Passenberg supplied their trajectories in ||6] for Fig. |8] through 
personal correspondence. 



TABLE IV 

1 1 -region autonomous switched system problem cost for 

Passenberg et al., traditional numerical programming 

(tnp), and tnp with switching point time iteration 

(TNP/ITER). 



Optimization 




Reported 


Cost 


Passenberg et al. 


m 


14.69 




Optimization 




Optimization 
Solution Cost 


Simulation 
Cost 


TNP 




14.71 


14.81 


TNP/iter 




14.69 


14.71 



the computed piecewise continuous controls of 14.81, 
slightly larger than the Passenberg et al. cost of 14.69. 

This solution can be improved by accounting for the 
(discontinuous) transitions across the switching mani- 
folds. To explicitly account for region transitions, we 
identified the time partitions in which they occur, es- 
timated the crossover times, appropriately subdivided 
each identified partition into two of unequal lengths, 
and then re-optimized. This resulted in a numerical 
programming cost of 14.69 and simulation cost of 14.71 
after 10 iterations on the estimated crossover times, 
which is consistent with Passenberg et al. reported cost. 
As an alternative, once the sequence of the discrete 
states is computed, methods that solve the hybrid optimal 
control problem assuming a known sequence of discrete 
states in, I25II , I31I can be used to compute the optimal 
switching points. 

TablellVllists the (similar) costs reported from Passen- 
berg et al. and those obtained here. The trajectory given 
by Passenberg et al. and the one obtained in the research 
reported herein using 172 uniform time partitions are 
plotted in Fig. [8j one observes the closeness of the 
trajectories computed by each of the two algorithms, 
Passenberg et al. and ours. The authors would like to 
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Fig. 8. Autonomous switched system witli 11 a-modes state path 
through the labeled state-space partitions: ( — ) simulation, (•) states 
from traditional numerical programming solution with 172 uniform 
time partitions, (' ') Passenberg et al. solution f6\. 



acknowledge that the algorithm of Passenberg et al. 
has an advantage due to the explicit handling of the 
switching- or discontinuity-manifolds. In order for clas- 
sical optimization such as SQP (as asserted in 15]) to 
obtain similar results, transitions across discontinuity 
manifolds must be properly accounted for Please note 
MPP, MIP/MPT, and CPLEX are not applied to this 
problem because they do not allow for switched PI 
penalty weights. 

III. Discussion 

The six hybrid optimal control problem examples 
presented demonstrate the differences and similarities 
between the embedding approach, MPP, MIP (MIP/MPT 
and CPLEX), and gradient based methods as well as 
the applicability of traditional numerical programming to 
autonomously switched systems. TablelV] lists key differ- 
ences between the various HOCP solution approaches. 
From the table, the embedding approach addresses a 
wider class of problems, does not require a specialized 
solver to implement, does not require additional variables 



for autonomous switches, and does not require supplying 
an initial mode sequence. Furthermore, the previous 
examples showed that the embedding approach led to a 
solution in all cases, the embedding approach solutions 
costs were less than or equal to the other methods costs 
(when available), and the embedding approach prob- 
lem was solved in less time (when available) than the 
other methods except for CPLEX. CPLEX's specialized 
routines resulted in solution times an average of one- 
third that of the embedding approach in the spring-mass 
example. The embedding approach utilizes MATLAB's 
finincon general purpose optimizer and is not streamlined 
for a specific class of problems Uke CPLEX. However, 
in the DC-DC boost converter example CPLEX took 
about 120 times longer than the embedding approach 
to solve the 12 partition control and prediction horizon 
boost converter problem, illustrating the "curse of dimen- 
sionality" associated with mixed-integer programming. 
The embedding approach solution can be implemented 
in real-time as illustrated in ||T3l where a boost converter 
problem is solvable in about 50 /is using a specialized 
algorithm. 

IV. Conclusion 

This paper delineates specific differences between the 
embedding approach for solving hybrid optimal control 
problems and multi-parametric programming, mixed- 
integer programming, and gradient-descent based meth- 
ods. First, the embedding approach's theoretical under- 
pinnings guarantee the existence of a solution under mild 
conditions. Second, the embedding approach generates 
a convex problem and is solvable with widely avail- 
able optimization packages such as sequential quadratic 
programming. Third, the embedding approach permits 
nonlinear state models as long as they are linear in 
the continuous time controls. Fourth, the embedding 
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TABLE V 

Hybrid optimal control problem solution approaches key differences (GDB is gradient descent based methods). 



Characteristic 



Embedded MPP MIP/MPT CPLEX GDB 0/0 



General nonlinear control models allowed 

Switched PI weights allowed 

Utilizes traditional numerical programming 

Computational complexity is NP-hard and rises exponentially with no. of modes 

Autonomous switches need discrete states assigned 

Problem initialized with predefined mode sequence 

Real time control Implementation 



Yes 


No 


No 


No 


YesA'es 


Yes 


No 


No 


No 


No/No 


Yes 


No 


No 


No 


No/No 


No 


Yes 


Yes 


Yes 


No/No 


No 


Yes 


Yes 


Yes 


-/- 


No 


No 


No 


No 


Yes/No 


Yes 


No 


No 


No 


No/No 



approach does not require any assumptions about the 
mode sequence. 



Appendix A 

Appendix: MATLAB Embedding Solution 

Approach Algorithm 



The embedding approach was applied to five hybrid 
control examples in the literature and results compared 
to the multi -parametric programming, mixed-integer pro- 
gramming, or gradient-descent method. The embedding 
approach produced equal to or lower performance index 
costs and faster online solution times than the other 
optimization methods with the exception of CPLEX 
in the spring-mass example. We attribute this to the 
specialized nature of the CPLEX solver as compared to 
the general MATLAB solver used to solve the embedded 
problem. Also, the embedding approach found control 
solutions when other methods failed. In the sixth hy- 
brid control example, an 11 autonomous-mode switched 
system which did not require embedding, we found the 
solution cost obtained with a hybrid minimum principle 
to be slightly less than that using traditional numerical 
programming; the important insight is that the numerical 
algorithm must account for discontinuities in the vector 
fields either implicitly or explicitly as is done in ||6l, Q. 
In summary, the embedding approach offers less problem 
set up and outstanding numerical results. 



Here we outline the steps to solve a hybrid optimal 
control problem with the embedding approach using 
MATLAB's /m/ncon function. 

1) (For continuous-time optimal control problem.) 
Transform the continuous-time embedded system 
dynamics into discrete time using the lengths of the 
prediction horizon partitions and either the forward- 
Euler, backward-Euler, or direct collocation with 
triangular basis functions ||9] methods. The result- 
ing discretized embedded optimal control problem 
dynamics are a series of nonlinear equalities. 

2) (For continuous-time optimal control problem.) Dis- 
cretize the PI cost over the prediction horizon 
partitions using trapezoidal numerical integration. 

3) Create equality constraints that enforce any 
supervisory-level interconnections and embedded 
mode value sum (embedded mode switch values 
sum to one) at the midpoint of each MPC partition 
if the problem was originally formulated in contin- 
uous time or partition boundary if the problem was 
given indiscrete time. 

4) Obtain an initial guess for the fmincon variables 
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to be solved for: states, algebraic variables, control 
inputs, and modes. The initial guess comes from 
the immediately previous EOCP solution, i.e., a 
warm start |32j, when t^ > <o- For the first EOCP 
solution, t^ ~ to, an in-house created prepro- 
cessing function is used to find the initial guess. 
This function works by first solving the EOCP (to 
a lower numerical tolerance than normal, llO^'^ 
versus 1-10^®) using fmincon for one partition 
ahead with a user specified initial guess. Then, 
partitions are successively added to the problem, 
with the EOCP solved after each addition, until 
the prediction horizon is reached. During partition 
addition, the solution associated with the previously 
added partition is carried forward to populate the 
initial guess of the current EOCP problem. 

5) Solve the EOCP using fmincon. 

6) Perform any mode and control projection on the 
controls for the first partition of the prediction 
horizon if the partition's EOCP solution is singular 
with some mode switch values in (0, 1). Theoret- 
ically, singular solutions occur when two of the 
possible mode-dependent Hamiltonians are numer- 
ically equal which occurs with small probability, 
generically, thus making bang-bang type solutions 
the norm. 

7) Either obtain a hardware measurement at the end 
of the current partition or execute a high accuracy 
simulation of the continuous-time switched system 
with the computed controls over one partition using 
MATLAB's ode23t function. 

8) Slide the starting time of the prediction horizon 
ahead one partition. If the end of the simulation 
is not reached, return to Step |4] 



Appendix B 

DC-DC Boost Converter Piece-wise 

Dynamics ifTTll 



The linear DC-DC boost converter dynamics for each 
duty cycle interval were developed using the procedure 
outlined in ifTTl . First, the duty cycle is divided into 
three regions: Dq = [0.0.45], L»i = [0.45,0.6], and 
D2 = [0.6,0.95]. Second, the allowable initial state 
conditions are specified. The scaled inductor current, 
i'l = ii/vs, initial condition is taken to lie in [0, 0.03] AA' 
while the scaled output voltage, v'^ = Vo/vg, initial 
condition is on [0, 1.6]. Third, the simulation data for a 
duty cycle interval is generated. The duty cycle, inductor 
current, and output voltage intervals are discretized using 
five evenly spaced points, then all possible combinations 
of the points are simulated over a time interval of tg 
duration using dZSl l. Mariehoz et al. |17'| do not provide 
the inductor current or output voltage intervals nor do 
they give the number of points in the discretization of the 
intervals. Fourth, the simulation data is fit to (l32l i using 
least squares. The resulting Unear models associated with 
the duty cycle intervals are 

0.9874 -0.01426 
0.2835 0.9939 



Ar. 



B„ 



Ar. 



B„ 



0.02828 
0.009680 



fn 



0.01347 
3.736-10-3 



(48) 



0.9862 -0.01184 
0.2342 0.9949 



B„ 



0.02000 
-0.009690 

0.9872 - 
0.1109 

0.01999 
-0.01526 



7 Jin,! — 

-5.600-10-3 
0.9965 



0.01431 
0.01074 



(49) 



Jn 



9.339 • 10-'* 
0.01477 



(50) 
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Appendix C 

Skid-Steered Vehicle Hybrid System 

Dynamics |fT9l 



Caldwell and Murphey |lT9l skid-steered vehicle 
mode 2 (front wheels sticking laterally and back wheels 
skidding laterally) and 3 (front wheels skidding laterally 
and back wheels sticking laterally) dynamic^ are given 
next. Mode 2 dynamics, f-z, are 



X 



F2 + F3 + Fi) cos(6'(t)) 



l-i-K-Mg . 

— sin(6l(t))v2,y 

-sm{e{t))[biFi-F2-F3 + Fi 



^j-KMga 



-V2. 



^ 12Jsm{9{t))[b{Fi -F2-F3 + F4) + a^^iKMg9[t) 

- aMX{t){nKgsm{e{t)) - cos(6'(t))^(i)) 
-aMY{t){pKgcos{0{t)) + sui{e{t))9{t))]/{al2Kj)\ 

- C2X{t) 



Y(t) = [X(t) +ta.n{e {t))Y{t)]e{t) 

+ ^7 sin(0(i)) tan{0it))[b{F, ^ F2 - F^ + F^) 
alvl 



a^KMg 



2 

tan(6l(t)) 
M 



«2,yJ 



[-{Fi+F2 + F3 + Fi)cos{0{t)) 



HKMg . 
^ — sin(6'(t))v2,j/] 

X sec(6i)[&(Fi -F2-F3 + F4) + a^fiKMge{t) 

- aMX{t){^Kgsm{e{t)) - cos{e{t))e{t)) 
+aMY{t){^iKg cos{e{t)) + sm{e{t))9{t))]\ 

- C2Y{t) 



(52) 



9{t) = -12[fe(Fi - ^2 - ^3 + ^4) + a^mMgm 
~ aMX{t){^iKgsm{e{t)) ~ cos{e {t))6{t)) 
+ aMY{t){^iKg cos{d{t)) + sin{e{t))9{t))]/{12Kj) 



with 



(53) 



(51) 



^Caldwell supplied the mode 2 and 3 dynamics through personal 
con'espondence. 



M =mb + Arrim 
J=4m,,„(a^+62) + 



Am,, 



12 



(wl + S-w^) 



"^(Bf+Bl) 



Kj=J + Ma^ 



V2.' 



sn\{e)X(t) + cos{0)Y{t) + ad{t). 



(54) 
(55) 

(56) 
(57) 
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Mode 3 dynamics, /s, are 

Xit) = jj[iFi+F2 + F3 + Fi) cos{e{t)) 
fJ-xMg . 

- i smie{t))[biFi - ^2 - F3 + F4) + ^^^^Y^V3,y] 
+ 12Jsm{0{t))[b{Fi - F2 - F3 + F4) + a'^fiKMgOit) 
+ aMX{t){nKgsm{e{t)) - cos(6'(t))^(i)) 
-aMY{t){fiKg cosie{t)) + siii{e{t))e{t))]/{al2Kj)\ 

- C3Xit) 



(58) 



Y(t) = [X(i) +tan(e'(i))r(i)]6'(i) 
1 



(aM) 
aUKMg 
2 

tan(6l(i)) 

a7 



sin(6l(t)) ta.n{e{t))[b{Fi - F2 - F3 + F4) 

[-(Fi +F2 + F3 + F4) cos(e'(t)) 



+ 2aA/12A-j {[^^^-^ + ^"'^^^ ~ 12^cos(20(t))] 
X scc(6'(i))[fe(Fi - F2 - F3 + F4) + a^fiKMge{t) 
+ aMX{t){^iK9sm{6{t)) ~ cos{e{t))6{t)) 
-aMY{t){jiKgcos{e{t)) + s[n{0{t))e{t))]\ 
- C3Y{t) 



(59) 



9{t) = -12[6(Fi -F2-F3 + Fi) + a^UKMgOit) 

+ aMX{t){^iKgsm{6{t)) - cos{e{t))e{t)) 

- aMY{t){pKg cos{e{t)) + sm{e{t))e{t))]/{l2Kj) 



(60) 



with 



Appendix D 

1 1 -Region Autonomous Switch Example 

Dynamics and State-Space Divisions |l6l, Q 

Passenberg et al. (|6l 11 -region autonomous switched 
system matrices {Aq, Bq), control bounds {Uq), PI ma- 
trices (5*^, Rq) and state-space partition line definitions 
are listed herq^- Autonomous-mode 1, Ui e [—2,2] x 
[-2,2]: 



(62) 





-1 


0.5 


^1 = 








-2 


-1 


51 = 


0.6 










0.6 





1 





B, = 











1 




1 





Wi = 











1 



Autonomous-mode 2, U2 £ [-2,2] x [-2,2] 



A2 = 


-1 -2 




0.5 -0.1 


82 = 


1 
1 





B2 = 


1 










1 


R7 = 


1 










1 



Autonomous-mode 3, Uz £ [—2,2] x [—2,2]: 



A, 



-0.5 -0.1 


-0.1 -0.1 


0.5 




0.5 





B3 



R3 



1 




1 




0.1 


0.1 



Autonomous-mode 4, U4 £ [-2,2] x [-2,2] 





-0.5 -3 


^4 = 






0.5 -1 




0.5 




^■4 = 


0.5 







1 




4 — 


1 






0.1 


4 = 


0.1 



Autonomous-mode 5, U5 £ [—2, 2] x [—2, 2] 



A. 



-2 
0.5 



B. 



1 
1 



(63) 



(64) 



(65) 



(66) 



(67) 



(68) 



(69) 



(70) 



V3^y = sm{e)X{t) - cos{9)Y{t) + a0(t). (61) 



^Passenberg supplied the problem values through personal corre- 
spondence. 
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0.1 
0.1 



i?5 = 



0.1 
0.1 



Autonomous-mode 6, Uq e [—2,2] x [—2,2]: 



Aa 



Se 



-0.1 0.2 


-2 -2 


8 




8 





Bfi 



Re 



1 








1 


1 








1 



Autonomous-mode 7, C/y G [-2,2] x [-2,2]: 





-0.5 


-2 


At = 








0.5 


-1 




1 




S7 = 


1 







1 




7 


1 






0.5 


7 


0.5 



Autonomous-mode 8, U^ £ [—1,1] x [—1,1]: 

-0.5 -1 
5 -0.5 



A, 



B. 



1 
1 



5*8 



0.1 
0.1 



i?8 



0.1 
0.1 



Autonomous-mode 9, Ug e [—2, 2] x [—2, 2] 



^c, 



59 



-0.1 0.5 


-0.5 -0.2 


10 




10 





Bg 



Rg 



1 








1 


1 








1 



(71) 



(72) 



(73) 



(74) 



(75) 



(76) 



(77) 



(78) 



(79) 



Autonomous-mode 10, Uw G [—10, 10] x [—1, 1] (u2 
dummy variable defined for consistency with other mode 
parameter definitions, has no effect on dynamics or cost): 



Sii — 



2 
2 



Rii — 



1 
1 



(83) 



Next, the lines separating the state-space regions are 
defined. 



xi = — 3, 


X2e[-4,-2] 


(84) 


xi = - 5, 


a-2e[-7,-3] 


(85) 


a:2 = - xi - 14, 


XI e [-7,-6] 


(86) 


2^2 = - 7, 


xi e [-10,5] 


(87) 


X2 — — 0.5x1 — 5.5, 


xiG [-10,-3] 


(88) 


a:2 = -4, 


•Tie [-3,5] 


(89) 


X2 =2x1 -4, 


XI = [-1.5,0] 


(90) 


a;2=0, 


XI e [-10,-5] 


(91) 


X2 ^ - Xi - 5, 


XI e[-5,-3] 


(92) 


X2=-2, 


.Tie [-3,0] 


(93) 


X2=-Xl-2, 


XX G [0,2] 


(94) 


X2 =a;i - 2, 


XI G [0,4] 


(95) 


2^2=2, 


xi G [4, 5] 


(96) 


X2 =0.8.Ti + 4, 


.TiG[-5,0] 


(97) 


a;2 = - 0.5x1 +4, 


XI G [0,4] 


(98) 



^1 



-4 


2 



2 


-4 


2 





-Bio 



i?i 



1 





1 





1 












(80) 



(81) 



Autonomous-mode 11, C/n G [-2,2] x [-2,2]: 



^11 = 



-1 1 
1 -2 



Bn = 



1 
1 



(82) 
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